  function [ysim,stsim,szero]=simssmodel(shocks,G,R,Z,A0,Pzero,C)
% function [ysim,stsim,szero]=simssmodel(shocks,G,R,Z,A0,Pzero,C)
% 
% SIMSSMODEL.M 
% 
% Simulate the path of the variables given the output of gensys.m 
% Inputs 
% ------
% shocks    [NX x T ] vector of shoclk
% G,R,Z     Output from gensys 
% A0        Mean of Initial State 
% Pzero     Variance of Initial State 
% C         Constant, Output from gensys
% 
% AJ 3/8/2005 Constant is in the observation, not in the state
% Initial state will be drawn 
% AJ 6/11/2009 szero is the initial drawn state ~N( A0, Pzero ) 
% ===================================================================
nz=size(Z,1);ny=size(G,1); 
if nargin == 6 || isempty(C)==1
    C=zeros(ny,1); 
end 
T=size(shocks,1);
shocks=shocks'; 
% Storage Matrices
ysim=zeros(nz,T);stsim=zeros(ny,T); 
C=Z*C; 
% Initial State 
%Pzero_chol=chol(Pzero); 
szero=mvnrnd(A0,Pzero); 
szero=szero(:); 
stsim(:,1)=G*szero+R*shocks(:,1);
ysim(:,1) =C+Z*stsim(:,1); 
% First Simulated obervation 
for ii=2:T ;     
    stsim(:,ii)=G*stsim(:,ii-1)+R*shocks(:,ii);
    ysim(:,ii)=C+Z*stsim(:,ii);  
end 
ysim=ysim';stsim=stsim'; 
